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Abstract 

An implicit method for the computation of 
unsteady flows on unstructured grids is 
presented. The resulting nonlinear system of 
equations is solved at each time step using an 
agglomeration multigrid procedure. The method 
allows for arbitrarily large time steps and is 
efficient in terms of computational effort and 
storage. Validation of the code using a one- 
equation turbulence model is performed for the 
well-known case of flow over a cylinder. A 
Detached Eddy Simulation model [1] is also 
implemented and its performance compared to 
the one equation Spalart-Allmaras Reynolds- 
Averaged Navier-Stokes (RANS) turbulence 
model [2]. Validation cases using DES and 
RANS include flow over a sphere and flow over 
a NACA 0012 wing including massive stall 
regimes. The project was driven by the ultimate 
goal of computing separated flows of 
aerodynamic interest, such as massive stall or 
flows over complex non-streamlined geometries. 

1. Introduction 

The use of time-dependent simulations for flows 
of practical interest is much less widespread than 
the use of steady-state flow simulations, due to 
limitations related to computational time and 
computational resources. However, many flows 
are inherently unsteady, particularly when large 
amounts of separation are present, and must be 
simulated as such in order to obtain meaningful 
results. 
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Transient flow simulations can either be based 
on explicit or implicit time-stepping schemes. 
Stability constraints limit the maximum 
allowable time step size of an explicit scheme, 
which is proportional to the smallest cell size in 
the entire computational mesh. Explicit schemes 
are well suited for unsteady applications in 
which the time scale of interest is comparable to 
the spatial scales. In these cases the mesh should 
be clustered only in regions of interest where it is 
absolutely necessary or the explicit time step can 
become unnecessarily small. 

Explicit time steps may become too restrictive 
for cases characterized by larger time scales or in 
cases where there is a wide variation in the grid 
resolution. In these cases it is desirable to 
develop a fully implicit method in which the 
time step is only determined by the physics of 
the flow and not by the cell size. At each time 
step the unsteady residual must be driven to zero 
and this is usually done using inner iterations. 
The number of inner iterations needed for each 
time step is related to the problem being solved 
and the size of the physical time step used. 

To develop an unsteady solver that can provide 
an accurate description in time of the flow is the 
first requirement to address before attempting to 
solve flows with massively separated regions 
that are inherently unsteady. On the other hand, 
it is also important to review the different 
techniques available to compute turbulence, 
which is at the core of separated flows, and will 
be a key parameter to determine how well the 
numerical scheme captures all the physical 
aspects of the flow. In this sense, numerical 
solutions of turbulent flow cases can be achieved 
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using different levels of approximation. The 
most widespread method is to solve the 
Reynolds averaged Navier-Stokes equations 
(RANS). In the RANS equations, the turbulent 
fluctuations appear in the Reynolds stress term 
that must be modeled using any of the turbulence 
models available in the literature. However, a 
common limitation of these models is their lack 
of generality since the model coefficients are 
usually set using simple well-documented flows. 
In this sense, current RANS solvers are fairly 
successful at predicting mostly attached flows, 
such as a wing in cruise condition, but fail to 
capture a range of different off-design situations 
as post-stall regimes, certain high-lift 
configurations, and non-streamlined bodies. In 
general, in the cases in which the RANS 
approach fails, the flow is characterized by large 
amounts of separation in which a very wide 
range of scales is present in the flow. While the 
small scales tend to be universal in nature, the 
larger scales are affected by the boundary 
conditions. This is the main cause of the lack of 
generality of turbulence models, as it is difficult 
to model the effect of the large scales in the same 
way for many different types of flows. 

The failure to develop a universally valid 
turbulence model has led to alternate approaches 
such as Direct Numerical Simulation (DNS) and 
Large Eddy Simulation (LES). DNS is the most 
straightforward approach to the problem. DNS 
consists in solving the governing equations on a 
mesh fine enough to capture the smallest scales 
contained in the flow with a scheme designed to 
minimize the numerical dispersion and 
dissipation. The drawback is extremely high cost 
of the DNS computation, which is proportional 
to at least Re 3 . For these reasons, DNS is 
generally limited to very simple flows and low 
Reynolds Numbers. 

The flow limitations of RANS and the difficulty 
of using DNS for realistic applied engineering 
problems have generated a great interest in the 
Large Eddy Simulation approach (LES) for 
computing flows with large amounts of 
separation. Large Eddy Simulation is a technique 
between Direct Numerical Simulation and 
Reynolds Averaged Navier-Stokes [3]. In LES 
the contribution of the large scales is computed 
exactly and only the smallest scales in the flow 
are modeled. However, a major difficulty of 
LES, particularly for external aerodynamic 
flows, is that near solid surfaces all the eddies 
are small and the “large” and “small” eddies tend 


to overlap. Therefore the required grid spacing 
and time step gradually fall towards DNS as the 
solid boundary is approached [4 ]. 

Using LES to resolve near wall streaks would 
bring an immense penalty at industrial Reynolds 
numbers. In this line of reasoning the Detached 
Eddy Simulation (DES) approach was conceived 
with the idea of combining the strengths of 
Reynolds Averaged methods near the solid 
boundaries and of Large Eddy Simulation 
elsewhere [1,4,5]. 

The ultimate goal of this work is to develop a 
large eddy simulation capability based on an 
existing unstructured grid Navier-Stokes solver 
[6] to be able to perform detached eddy 
simulations combing RANS near the walls and 
LES in massive separated regions in a non zonal 
manner and to compare the solutions obtained 
using this approach to the solutions obtained 
using a classical RANS approach with a Spalart- 
Allmaras one equation turbulence model [6]. 

2. Unsteady Reynolds Averaged Navier- 
Stokes Solver 

2.1 Steady solver description 

The Reynolds averaged Navier-Stokes equations 
are discretized by a finite volume technique on 
meshes of mixed elements, including tetrahedra, 
pyramids, prisms and hexahedra. In general, 
prismatic elements are used in the boundary 
layer regions, while tetrahedra are used 
elsewhere. 

Flow variables are stored at the grid vertices and 
a single unifying edge-based data-structure is 
used to handle all elements of the grid. 
Convective and viscous fluxes are discretized 
along edges using a central difference finite 
volume technique with added artificial 
dissipation. 

The non-dimensional steady conservative form 
of the full Navier-Stokes equations can be 
written as: 

R(w) = 0 (1) 

where w represents the solution vector of 
conserved variables and R represents the spatial 
discretization operator, or residual, which 
vanishes at the steady state. 
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An implicit solution procedure for solving 
equation (1) begins with the linearization of the 
residual about the current time step or iteration 
level: 


^ + R(w n )+^-Aw = 0 

3 r dw 


( 2 ) 


which can be solved as: 


The unsteady form of the governing equations is 
obtained by adding the time derivative to 
equation (1): 


— 0 (4) 

dt 

where R(w) denotes the discretization of the 
spatial derivative terms, as previously. 


•{-/?(»")} ( 3 ) 

where t is a pseudo-time used to advance the 
solution, and Aw represents the new correction to 
be applied to the solution. Rather than inverting 
the large Jacobian in equation (3), a reduced 
Jacobian which is simpler to invert is employed. 
In regions of isotropic grid cells, only the block 
diagonals of the Jacobian are retained, leading to 
a block-Jacobi or point-implicit scheme. In the 
boundary layer regions, where high grid 
stretching is present, lines are constructed in the 
direction normal to the solid wall boundary and 
the Jacobians along these lines are inverted, 
using a block tridiagonal solution algorithm. 
This procedure relieves the stiffness associated 
with high grid stretching in these regions, thus 
providing more rapid convergence [7 ]. 

This locally implicit scheme is used as the 
smoother on all levels of an agglomeration 
multigrid algorithm. The multigrid algorithm 
constructs coarse level grids automatically by 
fusing together neighboring fine grid control 
volumes to form a smaller number of large 
coarse grid control volumes. Convergence 
acceleration is achieved by cycling back and 
forth between the fine and coarse grid levels of 
the multigrid sequence, using the locally implicit 
solution procedure as a solver on each grid level 

[7L 

The unstructured multigrid solver is parallelized 
by partitioning the domain using a standard 
graph partitioner [8] and communicating 
between the various grid partitions running on 
individual processors using the MPI message- 
passing library [9]. The solver can be run on 
distributed or shared memory machines, 
including clusters of personal computers. 

22 Unsteady term implementation 


Making use of a second-order accurate three 
point backward approximation for the time 
derivative and evaluating R(w) at time level 
(n+1), we obtain [10,11,12]: 

— -wT 1 .W+ — . w /'- I +/?( H ' ,+1 ) = 0 (5) 

2 At At 2 At 


We may now define and unsteady residual as: 

R* (w) — — + R(w) — 5 (h |B , H’ n_l ) (6) 

2- At 

with 

S(w",w n ~') = — w n — w nl (7) 

At 2- At 

where w/ l+l is the approximation to w and the 
source term S(w n ,w /l1 ) remains fixed throughout 
the solution procedure at each time step. In this 
form, the non-linear problem to be solved at each 
time-step reduces to R*(w) = 0. Using this 
redefinition of the residual, the same iterative 
multigrid procedure employed to solve the 
steady-state problem can be used to solve the 
non-linear unsteady residual at each physical 
time step [10,1 1,12]. 

The solution procedure consists of an outer loop 
over the physical time steps, used to advance the 
problem in time, and an inner (multigrid) loop 
over pseudo-time used to drive the unsteady 
residual to zero at each time level. While the 
maximum size of the physical time step is only 
determined by the physics of the problem, the 
convergence of the inner sub-iterations for each 
time step varies inversely with the size of the 
physical time step. 

3. Unsteady solver validation 

The flow around a circular cylinder is a well- 
known case, which has been widely studied 
computationally and experimentally. This case is 
used as the basis for validation of the unsteady 
RANS solver, and for assessing grid resolution 


A >v = 


A r dw 
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and time step requirements for accurately 
predicting the vortex shedding frequency 
observed in the cylinder flow. Two different 
meshes of 252,000 and 631,000 grid points and 
three different time steps of 0.5, 0.25 and 0.1 
were used. The time is non-dimensionalized as t 
= to/(d/Uoo) where d is diameter of the cylinder 
and Uoo is the ffeestream velocity. 



Figure 1: 3D View of cylinder and lateral walls 

The one equation Spalart-Allmaras turbulence 
model [2] was used for all calculations in fully 
turbulent mode. In all cases the agglomeration 
multigrid strategy was used with four levels. The 
Mach number is 0.2 and the Reynolds number is 
1200 for this case. All runs were performed in 
parallel using 16 processors of a Pentium PC 
cluster at ICASE. 

The computational domain in the plane normal 
to the cylinder span has an aspect ratio of l and a 
side length of 100 cylinder diameters. A span of 
two cylinder diameters is employed, and inviscid 
(slip velocity) boundary conditions are applied at 
the end-walls. The three-dimensional 
simulations reported herein were also compared 
with two-dimensional simulations of flow 
around a circle using a validated two- 
dimensional unstructured solver [12], and found 
to agree well in terms of force coefficient 
histories and shedding frequency. 

Table 1 shows the Strouhal Numbers computed 
for each mesh and each time step of the three- 
dimensional simulations. Convergence is 
achieved as the time step is reduced and the 
mesh size increased. A second-order accuracte 
convergence behavior is observed as the time- 
step is reduced, validating the accuracy of the 


three-point backwards difference scheme used to 
discretize the time step. From the smallest time 
step results, the solution can be seen to be grid 
converged, at least with respect to the prediction 
of the vortex shedding frequency. The computed 
Strouhal number compares very well to the 
experimental value of St = 0.21 given by Dresher 
[13,14]. Figure 2 shows the time history of the 
lift coefficient, while the oscillatory pattern 
corresponding to the vortex shedding is shown in 
Figure 3. 



Time Step 

0.5 

0.25 

0.1 

% 

<3 

0.252 

Million 

0.19249 

0.20304 

0.20833 

0.631 

Million 

0.19379 

0.20408 

0.20833 


Table 1, Predicted Strouhal Number for Various 
Grid and Time Step Sizes 



Figure 2: Mach contours at two different time 
states for flow over circular cylinder . Mach = 
0.2, Re = 1200. 
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Figure 3 . Computed LiftCoefficient Time History 
for Flow over Circular Cylinder using three 
different time steps. Mach = 0.2 , Re = 1200 


Mean pressure coefficient distribution over the cylinder. 
M * 0.2, Re * 1200 



Figure 4 . Mean Surface Pressure Coefficient 
Distribution for Flow over Circular Cylinder at 
Mach = 0.2 , Re = 1200 , Compared with 
Experimental Data at Various Reynolds 
Numbers 


Figure 4 shows the pressure distribution over the 
surface of the cylinder. The pressure distribution 
was computed by averaging results at different 
times distributed along several oscillations. For 
this calculation the finest mesh and the smallest 
time step was used. As can be observed from 
Figure 4, the computed pressure distribution 
compares closely to experimental results at a 
higher Reynolds number than the one used for 
these computations. This is likely due to the use 
of the turbulence model in fully turbulent mode, 
in order to avoid the issues of transition 
prediction. Similarly, the backpressure obtained 
of C pb = -1.20 compares closely to the C pb = - 
1.20 at Re=27,700 measured by Linke [13,15]. 
The mean value of the drag coefficient computed 
is C d =1.3 compared to C d =1.2 as measured by 
Wieselsberger [13,16] for Re=30,000. 


4. Detached Eddv Simulation. 

Detached Eddy Simulation (DES) is a hybrid 
technique that combines RANS and LES in a 
non-zonal manner. DES is based on the Spalart- 
Allmaras one equation RANS model [2] in 
which the length scale d, which is traditionally 
taken as the nearest distance at any given point to 
the closest wall, is replaced as the minimum 
between the distance to the wall and a length 
proportional to the local grid spacing: 

d DES = min (d, C DES Ax) 

where C DES represents a model constant which 
has been taken as 0.65 in previous work [1,17]. 
Traditionally, on structured grids [1,17,18], Ax is 
taken as the maximum grid spacing over all three 
directions. In our particular case, the definition 
of Ax has been modified for unstructured grids 
by taking it as the maximum edge length 
touching a given vertex. In boundary layer 
regions, Ax far exceeds the distance to the wall 
d, and the standard Spalart-Allmaras RANS 
turbulence model is recovered. Flowever, away 
from the boundaries d exceeds C DE s Ax and the 
models turns into a simple one equation sub- 
grid-scale (SGS) model with the mixing length 
proportional to the grid spacing. This effect is 
illustrated by plotting contours of the distance or 
length scale function for both the RANS and 
DES models in Figure 5, where it is observed 
that both models employ the same length scales 
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near the wall, but use vastly differing length 
scales in the regions far removed from the wall, 
where the DES model reverts to an LES mode 
and a Smagorisnky-like expression for the eddy 
viscosity is obtained. 


URANS DES 



Figure 5. Distance function/Length Scale 
comparison between URANS and DES . 


5. Validation. 

5.1 Flow around a Sphere. 

DES is applied to predict the flow around a 
sphere. Similarly to the flow around a circular 
cylinder, this case has been widely studied using 
experimental and computational approaches. It 
has been shown that over a wide range of 
Reynolds numbers (280< Re < 3.7x1 0 5 ) the flow 
is characterized by the vortex shedding with 
large-scale vorticity emanating from the shear 
layer which separates from the surface of the 
sphere. In our case a Reynolds number of 10 4 is 
used, corresponding to the sub-critical regime 
(laminar boundary layer separation). This is 
similar to the computations performed by 


Constantinescu [17], which provides a 
comparative basis for our results. 

The main goal of this case is to compare the 
results obtained using traditional URANS and 
DES. The comparison includes the time history 
of integral parameters such as drag coefficient 
and mean distribution of pressure around the 
sphere. 

5.2.1 Numerical approach. 

The flow over a sphere is calculated for a 
Reynolds number of 10 4 and a Mach number of 
0.2 on an unstructured mesh of 767,000 vertices. 
The computational domain is a cube of 100 
sphere diameter lengths in each direction. This 
large domain is chosen to ensure that the 
downstream boundary condition does not lead to 
spurious oscillations within the domain. Near the 
surface of the sphere, the normal grid spacing is 
l O ' 4 x d, where d represents the diameter of the 
sphere. 

For the unsteady calculations two different time 
steps of 0.1 x d/Upo and 0.05 x d/U«, were used. 
In all cases the multigrid strategy was used with 
four grid levels. All runs were performed on 32 
processors of a Pentium PC cluster at ICASE. 

For the DES runs, the C D es value used was 0.65. 
This value was chosen based on previous studies 
by Shur et al. [1], where the value was calibrated 
in isotropic turbulence. Moreover, 
Constantinescu et al. [17] prove this value of 
C D es as the most appropriate for their 
calculations of the flow over a sphere at Re = 
10 4 . A parallel study is being conducted to 
examine the effects of variations in the value of 
C D hs using the same test case of decaying 
homogeneous turbulence in a box as done by 
Shuretal. [1]. 


5.2.2 Code validation. 

Following the procedure used by Constantinescu 
et al. [17], the accuracy of the numerical 
approach was established by comparing the 
results obtained in the steady regime at Reynolds 
number 250 with previous computational and 
experimental results. 

The drag coefficient is compared in Table 2 to 
other simulation results and experimental data. 
Because previous results were based on 
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incompressible simulations, and the current 
solver is a density based compressible 
formulation, the importance of compressible 
effects was also investigated by running the 
simulation at Mach numbers of 0.2 and 0.1, both 
with and without a low Mach number 
preconditioner [18,19]. The agreement is 
satisfactory for all the cases tested. 



Cd 

Constantinescu et 
al. 

(2000) 

0.70 

Johnson and Patel 
(1999) 

0.70 

Experimental 

0.70-0.72 



M = 0.1 

0.7141 

M = 0.2 

0.7015 

M = 0.1 
Low Mach 
Number pre- 
conditioner. 

0.6961 

M = 0.2 
Low Mach 
Number pre- 
conditioner. 

0.6950 


Table 2. Computed Steady Drag Coefficient for 
Flow over Sphere at Re = 250 compared with 
Experimental and Previous Computational 
Values 


5.2.3 Unsteady Results. 

For the unsteady runs, the flow around a sphere 
is computed at a Mach number of 0.2, without 
any additional low Mach number 
preconditioning, and a Reynolds number of 10 4 , 
At this Reynolds number the detached vortex 
sheet from the sphere is fully turbulent while the 
boundary layer on the sphere remains laminar. 
The Strouhal number associated with the vortex 
shedding at this Reynolds number is in the range 
of 0.185-0.200 depending on the investigation. 
The large scatter of the data is mainly due to the 
influencing parameters and the measuring 
methods of the different investigations. 



Time 


Figure 6 . DES Time history of Drag Coefficient 
for flow over Sphere at Mach = 0.2, Re = 10,000 
using two different time steps 



Figure 7. URANS Time history ’ of the drag 
coefficient for flow over Sphere at Mach = 0.2 , 
Re = 10,000 using two different time steps 
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Figure 8. Comparison of URANS and DES 
Mach Contours for flow over sphere. 


Moan pressure distribution over the sphere. 
M- 0.2. Re -10.000 



Figure 9 .. URANS and DES Average Surface 
Pressure Distributions Compared with 
Experimental Data at Various Reynolds 
Numbers for flow over Sphere 
Mach = 0.2, Re = 10,000s 

The time history of the drag coefficient is shown 
in Figures 6 and 7 and reveals important 
differences between URANS and DES. The 
mean value of the drag coefficient in both cases 
is close to the experimentally reported value of 
0.40 from Schlicting [20]. However, the 
frequency content in each case is completely 
different. The URANS simulation appears to 
damp out most of the oscillations present in the 
DES run, while the DES runs show a very 


chaotic oscillatory pattern quite similar to the 
solutions obtained by Constantinescu et al [17]. 
Spectral analysis of the time -dependent drag 
coefficient history reveals a peak corresponding 
to a Strouhal number of 0.1 which is not in 
agreement with the values 0.18 to 0,2 reported 
experimentally [21]. This may be due to an 
insufficiently long time history sample, since 
less than three full periods of this frequency are 
present in our sample. 

The Mach number contours depicted in Figure 8 
corroborate the difference in the predicted flow 
using regular Spalart-Allmaras RANS turbulence 
model and detached eddy simulation (DES). 
DES shows a wider range of scales present in the 
flow while regular RANS models tend to smooth 
out the smaller scales. This is the effect expected 
from DES since the length scale redefinition 
increases the magnitude of the destruction term 
in the Spalart-Allmaras model, drawing down 
the eddy viscosity and allowing instabilities to 
develop. Predictions of the mean pressure 
distribution over the surface of the sphere are 
shown in Figure 9. The surface pressure is seen 
to match well with experimental results at 
Re= 165,000 reported by Achenbach [21], which 
is in agreement with the results reported by 
Constantinescu [17]. The average computed 
separation angle of 81° compares resonably to 
experimental vale of 82.5°. 

5.3 Flow over a NACA0012 airfoil. 

RANS and DES were used to compute the flow 
over a NACA0012 airfoil at a Reynolds number 
of 10 5 and a Mach number of 025. Runs were 
performed for 10 different angles of attack 
ranging from 0 to 16. Our main interest was to 
study and compare RANS and DES in the stall 
regime. Other studies (Shur et al. [1]) have 
focused their interest in higher angles of attack, 
such as 45, 60 or even 90 degrees, showing very 
good and promising results in which DES is 
within 10% of the experimental values, while 
traditional RANS predict results with 40% 
errors. All runs are computed as fully turbulent 
to avoid having to trigger the transition point. A 
mesh of 966,00 node points is used with a first 
grid spacing normal to the solid boundary of 10 5 
chords. The spanwise boundary conditions are 
slip (inviscid) walls, and the span of the 
geometry is two chords lengths. 

In this case, due to time and computing 
constraints, only one time step of 0.25 x cAL, 
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was used, where c is the airfoil chord. The 
calculations were carried out for a maximum of 
50 time units for the post-stall cases, where a 
time unit represents the time it takes for the 
undisturbed far-field flow to travel one chord 
length. 

5.3.1 Results. 

The time history of the lift and drag coefficient 
show good agreement between URANS and 
DES for angles of attack bellow 11.5 degrees. 
This was expected since for pre-stall conditions 
the DES model operates primarily in the RANS 
mode. However, for the post-stall condition, i.e. 
angles of attack over 1 1.5, the time history of the 
force coefficients obtained using URANS and 
DES showed differences similar to the ones 
observed for the sphere case. The time history of 
the DES results shows higher unsteadiness that 
URANS, indicating that more scales are being 
captured in the separated region. The DES 
simulation predicts a more severe stall than the 
URANS results, i.e. lower post-stall lift and 
higher drag. However, both computations stall at 
the same incidence and comparison with 
experimental data [22] is not sufficiently close to 
favor agreement with one method over the other. 
This is in contrast to previous work [1] where 
DES was reported to provide closer agreement 
with experiement for very high angles of attack, 
such as 45 to 90 degrees. However, little 
attention was devoted to the prediction of the 
onset of stall in [1]. 




Time 


Figure IL Sample DES Drag Coefficient Time 
History for Flow Over NACA00J2 Wing at 
Various Angles of Attack 


Urtooefflctent NACA0012 
R«=1O J ,M = 0JZ5 



Figure 12 . Comparison of Computed Lift 

Coefficient versus Angle of Attack for URANS 
and DES versus Experimental Data [22] at two 
Different Reynolds Numbers 


Figure 10: Sample URANS Drag Coefficient 
Time History for Flow over NACA00I2 Wing at 
Various Angles of Attack 
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Drag coefficient NACA0012 
Re - 10\ M = 0,25 



Figure 13.. Comparison of Computed Drag 
Coefficient versus Angle of Attack for URANS 
and DES versus Experimental Data [22] at two 
Different Reynolds Numbers 


6.Conclusions. 


A steady-state unstructured multigrid solver was 
extended to unsteady flows and validated by 
computing the flow over a circular cylinder. 
Results show good agreement with experimental 
results in Strouhal number, pressure distribution , 
and flow patterns over the cylinder. 

A Detached Eddy Simulation (DES) technique 
was implemented in the code and compared with 
regular URANS by computing the flow over a 
sphere and the flow over a NACA0012 airfoil. 
The results for the sphere show that DES 
provides a better description of the flow in the 
separated regions. This was concluded from the 
time history of the force coeffcients and from the 
pressure distribution over the surface of the 
sphere. The time history of the drag coefficient 
computed using DES showed small scales 
present in the flow that URANS did not capture. 
The Strouhal number obtained does not match 
the experimental value but longer time histories 
are required for a more definite determination of 
this values. 


Prediction of the onset of stall for a NACA0012 
wing was not substantially improved using the 
DES approach over the URANS approach. 
However, DES has been reported to provide 
better agreement with experiement at very high 
angles of attack, and this will be investigated in 
future work. Additional work is also underway 
to simulate the decay of homogeneous isotropic 
turbulence in a periodic domain, in order to 
calibrate the DES model for unstructured grids. 
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